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An extended Kalman filter is used to estimate the 
velocity of an object moving across an image frame and to 
reduce the undesirable effects of noise. The extended Kalman 
filter is implemented in the spatial frequency domain to 
reduce the number of computations. The resulting filter 
structure is a parallel bank of third-order extended Kalman 
filters. This parallel structure is referred to as the 
modified extended Kalman filter. The performance of the 
modified extended Kalman filter is evaluated under a variety 
of noise conditions using computer simulations. Simulations 
employed two test objects moving across a reference image in 
the presence of zero-mean, white, Gaussian noise. The 
performance of the filter was demonstrated when these objects 
were moved at integer and noninteger velocities. Performance 
was also evaluated when a stationary background was included 
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THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in 
this research may not have been exercised for all cases of 
interest. While every effort has been made within the time 
available to ensure that the programs are free of computational and 
logic errors, they cannot be considered validated. Any application 
of these programs without additional verification is at the risk 


of the user. 
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I. INTRODUCTION 


In many radar and remote sensing operations, successive 
image frames may contain a moving object of interest. 
Identification of this object may be complicated by image 
degradation resulting from the presence of background noise. 
The idea of using an extended Kalman filter (EKF) to enhance 
the quality of the image frame and provide an estimate of the 
object velocity was proposed by Burl (Ref. 1]. 

In this thesis, a moving image model is defined as a 
series of image frames containing a moving object. An image 
frame consists of a two-dimensional array where the individual 
array elements are defined as pixels. Each pixel is given a 
numerical value based on the intensity of the image at that 
point. These values are either stored in a computer or 
projected on a video display terminal. 

The EKF requires a model that describes the evolution of 
the image frames in time. This spatial domain model is 
derived in Chapter II. A shift operator in the model 
describes the motion of the object as it traverses the image 
frame and is a function of the object’s velocity. A 
prohibitive number of calculations would be necessary if the 
EKF were implemented in the spatial domain. [Ref. 1] 

Chapter II further describes the moving image model after 


it has been transformed to the spatial frequency domain. The 


two-dimensional discrete Fourier transform (DFT) transforms 
these two-dimensional image frames from the spatial domain to 
the spatial frequency domain. By transforming the image frame 
to the spatial frequency domain, the EKF can analyze each of 
the spatial frequencies separately. The moving object is 
evidenced by a phase shift in each of the spatial frequencies. 
The implementation of the EKF in the spatial frequency domain 
reduces the number of calculations. The number of 
calculations can be further reduced by llmiting the number of 
spatial frequencies that are analyzed. [Ref. 1] 

Chapter III presents the extended Kalman filter equations. 
As the velocity estimate error approaches zero, the EKF is 
shown to converge to a parallel set of third-order extended 
Kalman filters. This parallel structure is referred to as the 
modified extended Kalman filter (MEKF). (Ref. 1] 

Chapter IV discusses the results of the computer 
simulations that were run to record the performance of the 
EKF under a variety of noise conditions and background 


environments. Conclusions are presented in Chapter V. 


II. THE MOVING IMAGE MODEL 


A. THE MOVING IMAGE MODEL IN THE SPATIAL DOMAIN 
Successive image frames containing a moving object can be 
modeled by a nonlinear state equation in the spatial domain 


as, 


E. | ۶)1+1( | | 8(v) 0 | | N, | (kK) 


v(k+1) 05 v(k) t (K) 


= a(v)x(k) + ۵ )16( . 1) 


An image frame has the dimensions of N; by N, pixels. The 
amplitude of each pixel, when stacked into a column, forms 


the vector f, where f e RU, 


The velocity vector, v, 
consists of a two-element column vector describing the 
horizontal and vertical velocities of the moving object in 
terms of the number of pixels per iteration. The vectors f 
and v combine to form the state vector x, where x c ې‎ 129 
S(V) is a two-dimensional shift operator with the magnitude 
and direction being a function of the velocity vector, v. 
The plant noise vector, ¢, is composed of the image plant 
noise, Ç r and the velocity plant noise, ¢,. The plant noise 


1s assumed to be a zero-mean, white, Gaussian random process 


With a covariance matrix defined as 


o I 9 
[ 5, 5, ] : = > ; (2.2) 
DEC 


where E[e] is the expectation operator and I is the identity 








operator. 
The measured image, y(k), is defined by the corresponding 


measurement equation, 


y(k) = f(k) + v(k) 
= { I0 ] x(k) + v(k) 


= ec x(k) + v(k) , )27 


C2) (nana) 


where y € R and v € R The zero mean, white, Gaussian 


random process, v(k), has a known covariance matrix 


R = E [ u(k) n! (K RIS po IM (298p 


VN U 


B. DERIVATION OF THE SHIFT OPERATOR IN THE SPATIAL FREQUENCY 
DOMAIN 


A image frame is a two-dimensional, finite-extent sequence 
defined as f(ni,n)) where O < rn; < Ni-1 and O < n < N51. The 
sequence f(ni,n,)) is transformed from the spatial domain to the 
Spatial frequency domain using the two-dimensional discrete 


Fourier transform (DFT), 





E ۱ 1 Dos 
F(K ,K,) = > 2> f(n ,n)exp لرل وغل‎ , (259 
N, N, 


n, -0 1 - 0 


where 0 < ۳ > 1۱-1 and 0 > <#N2-1. 

The shift operator in the spatial frequency domain is 
derived from the circular shift property of the discrete 
Fourier transform [Ref. 2:p. 67]. The circular shift property 


is defined as 


mk 


x( (O7m)), ,(CGo7m)), ) e Wy, وو‎ X (X, (ئ1,‎ 


e D(m,,m)X(k,,k;), (2.6) 
Ec 
N 
where تو )اا‎ is the modulus operator, Wy, = e Í and D(m,,m,) 
i 
1s the transformed shift operator. Equation 2.6 indicates 


that a circular shift in the spatial domain results in a phase 
shift in the spatial frequency domain. 

The original image frame, x(n,,n), is assumed to be a 
periodic two-dimensional sequence with the fundamental period 
equal to the frame dimensions. The object moving across an 
image frame then describes a circular shift. [Ref. 2:pp. 61- 
62]. This implies that, as the object moves off one edge of 
the screen, it reenters the screen from the opposite side. 
This assumption was necessary because the processing 
capabilities of the computer used to execute the EKF algorithm 
limited the image frame dimensions. Image frames of finite 
extent, such as a radar screen or video display, do not 
require this property because the moving object’s size will 
likely be  disproportionally smaller than the screen 


dimensions. 


Derivation of D(m,,m) begins by introducing a circular 
shift into the original image frame, f(n,,nj). Equation 2.5 


now becomes 


21۳۱ 2Tn5k»5 
N,-1 N5-1 EU "No NES 
ni = 0, n. = Ü 


where F'(k,,k?) is the circular-shifted, transformed image 
frame. The amount of pixel movement is described by m, for 


the horizontal shift and m; for the vertical shift. 


Letting 1, - n,-mj| and 1; - n;-m;, Equation 2.7 becomes 
2m(1,* m,)k, —2mT(15* m;)k 
Nic: “AlN, = ملا کا سای‎ 
BA (k k) =" > f(1,,1,) e . (209 
تزا‎ 0 o c 


Rearranging terms, Equation 2.8 reduces to 


5 2Tm k, ۲ 2Tm^k» 
l N N2 


FP f (k K.) mm m PORE 
= ا ا‎ 2 F(X,, X;) 
— D(m,,m) F (Ki; K) , (2.9) 


where F(k,,k2) is the transformed image frame without the 
circular shift. This proves that, when an object moves across 
an image frame, the motion is described by a phase shift in 
the spatial frequency domain, thereby agreeing with the 


clrcular shift property (Eqn. 2.6). 


The notation is simplified by defining the frequency 


vector, o, as 





W, N, 
— سم‎ 2 ۰ 10 
۴ دك‎ 2nk, 


and the velocity vector as 


۷ m 
v-a J-la]. C 
V- m, 


Equation 2.9 now becomes 


ET 
HE راع لی‎ 


D(m,,m,) F(K,,K,) 


Diy) ETC ٦ یج‎ 


where T is the transpose operation. 
The example contained in Appendix A further demonstrates 
the use of the phase shift operator in reconstructing image 


frames in the spatial frequency domain. 


C. THE MOVING IMAGE MODEL IN THE SPATIAL FREQUENCY DOMAIN 
The spatial domain model described by Equations 2.1 and 

2.3 can now be transformed to the spatial frequency domain. 

The transformation is accomplished by using the two- 


dimensional discrete Fourier transforn, 


E kK) = ><>, > f(n,,n5)ex -j2 8 are 29] 3 
| هئ‎ EN i P Jen N un N. , ( ° ) 


n, -0 ۲ -0 


where 0 > k, < N4-1 and 0 > ۶ «x N;-1. 


The vector, f(k), and its corresponding plant noise, سو‎ 


are transformed to the spatial frequency domain, 


Re (F, (kK) ) 
Im(F,(k)) 


F(k) » F(f(k)) - ME (2.14) 


1 Fu iw, (X) ] 
2 
v (k) 


Z(k) = F(t,(kK)) , (2309 
where F(e) is the two-dimensional discrete Fourier transform. 
(It should be noted that the index k refers to the sample at 
time k and is not related to the frequency vector.) Only half 
of the spatial frequencies need to stored in F(k) due to the 
conjugate symmetry that results from the transformation of a 
real image (Ref. 3:p. 45]. 

The state equation is defined in the spatial frequency 


domain as, 


X(k+1) = A(v)X(k) + 2(k) 
F(k+1) D(v) 1 MN Zp (Kk) 
وس‎ p v (X) ال‎ — 


(n4n2*2) 


where X € R The transformed shift operator, D(v), is a 


diagonal matrix where the diagonal terms, D;(v), are the 


transformed shift operators for a specific spatial frequency. 
Because the real and imaginary parts of the spatial frequency 
vector are stacked into a column vector, D;(v) is expressed 


in terms of its trigonometric relations, 


F,(k+1) = D,(v)F, (k) 


Re (F, (K+1) ) | cos(w,'v)  sin(w,'v) | Re (F, (k) ) 
= ; T T 5 (2:17) 
Im(F,(X*1)) -sin(@, V) cos(o, V) INE (KY) 


(N -1) (N.-1) 
2 


where i = 0,1,..., 


The covariance matrix for the image plant noise, Z, is, 


6071 0 
پ‎ 2 ۶ ۱ i | ۱ 0.187 


O oI 


۷ 


where E[e] is the expectation operator and I is the identity 
operator. Because the discrete Fourier transform is a linear 
operation on dm Z remains a zero-mean, Gaussian random 
process. 

The corresponding measurement equation in the spatial 


frequency domain becomes, 


Y(k) = F(k) + N(k) 


[ IO ] X(K) + N(K) 


C X(k) + N(k) (2.19) 


. > 


n E 


where Y(k) = F(y(k)) and N(k) » F(v(Xk)). The state equations 
given by Equations 2.16 and 2.19 are the basis for 


implementing the EKF in the spatial frequency domain. 
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III. THE EXTENDED KALMAN FILTER 


A. GENERAL 

In this chapter, the extended Kalman filter (EKF) is 
obtained from the state equations that were derived for the 
moving image model. Equations 2.16 and 2.19 are used to 
implement the EKF in the spatial frequency domain. The 
modified extended Kalman filter (MEKF) is based on the 
diagonal properties of the transformed shift operator, D(v), 
and covariance matrices. The parallel structure of the MEKF 
is formed in the limit as the velocity estimate approaches 


the actual velocity. 


B. THE EXTENDED KALMAN FILTER 

The extended Kalman filter (EKF) algorithm is used when 
the state equations are linearized about each new estimate as 
they become available [Ref. 4: p.158 ]. The EKF uses the 
spatial frequency domain state equations, Eqns 2.16 and 2.19, 
for the moving image model. 


The linearized state equation is given as 


6X(k+1) = A(X,)8X(k) + Z(k) , (Soll) 


where §6X(k)=X(k)-X,. Here, X, is shorthand notation for the 
current estimate X(k|k) [Ref. 4: p.158] and A(X,) is the 


Jacobian of the nonlinear state equations about X,. This 
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means that the ij component of Á(X,) is the partial derivative 
with respect to the jth state of the ith component of the state 


equation (Eqn. 2.16) and Á(X;) is then composed as follows: 


3 N.N 
A(X,) = diagonal([ (X. )] 1 K (3.25 





cos(®) sin(a) [-f,(k)sin(a«)+f,(k)cos(«) Jo’, 


A, (X,)= -sin(a@) cos(«) [-f,(k)cos(«)-f,(k)sin(a) Jw’, |» (3-3) 
0 0 1 
where e = o'v, f (k) = Re[F,(k)] and £,(k) = Im[F, (k)]-. A 


further discussion of the linearization techniques used by 
EKFs is contained in Ref. 4, p.154, and Ref. 5, p.194. 
The EKF equations for the moving image model are 


summarized below. 


Prediction 
State prediction 
X(kTI|K) = AV) ae) (3.4) 


Covariance prediction 


P(k+1|k) = A(X,) K(k] k) P(kIK)A (X,) X(k] x) به ې‎ )3 


Innovation 
e(k+1) = Y(k+1) ~ YV(Kk+1|K) 
= ¥(k+1) - C X(Kk*1|k) (3.6) 
Kalman gain 
G(k*1) 2 P(k«1|k)C'[C P(kri|]|k) C! & RJ! (3.2) 
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Correction 
State correction 
X(k+1]k+1) = XK(k+1|k)+G(k+1) [Y(k+1) -Y(k+1|k) ] (3.8) 
Covariance correction 
P(k+1|k+1) = [ I - G(k+1)C]P(k+1|K) (3.9) 
The notation 1۳11۴ reads as "at time k-*1, given data up 


through time k." 


B. MODELING PROBLEMS 

When the Kalman filter is implemented, the operation of 
the filter can be degraded by modeling errors resulting from 
linearization around the imperfectly estimated state. One 
possible method of reducing this divergence is to add 
fictitious plant noise to the system model. (Ref. 6: p.280] 

The system dynamics matrix, A(v,+A), is assumed to have 
model errors which is represented by A in Eqn. 3.10 and is 
expressed as a first-order Taylor series expansion. The 
first-order term is then combined with the actual plant noise 


vector, w(k), to form the plant noise vector, ¢. 


X(k+1) = A(v,t+4)X(k) + w(k) 


= A(v)X(k) + ۱ I "PE X(k) + w(k) | 





= A(v)X(k) + t(X) , (3.10) 


where A denotes the model error and w(k) is the actual plant 


noise vector. Including the model error, as part of the plant 
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noise vector, prevents new measurements from being weighted 
too lightly because the Kalman gains tend to zero as t oo. 
This ensures that the model continues to track the system as 


new measurements are taken. 


C. THE MODIFIED EXTENDED KALMAN FILTER 

As mentioned earlier, the diagonal properties of the 
transformed shift operator, D(v), and the covariance matrices 
are the basis for the parallel structure shown in Figure 3.1. 

This parallel structure is referred to as the modified 
extended Kalman filter. Reference 1 outlines how the EKF 
converges to the MEKF as the velocity estimate approaches the 
actual velocity. 

Practical implementation of the MEKF suggests that each 
filter in the parallel bank be dedicated to a specific spatial 
frequency. Each individual filter is referred to as a single 
frequency extended Kalman filter (SFEKF). 

The state vector for a specific spatial frequency is 


defined as: 


X S) Re(F (K) ) 
X (k) = | X, (K) = | Im(F,(k)) , (3.11) 
X, , (K) w,'v 


where the first two states, X;,(k) and X; 2(k), are the Four ٢ ٢ 
coefficients associated with a specific spatial frequency and 


the third state, X;3(kK), is the dot product of the spam 
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frequency vector and the velocity vector. The reason ۰ 
making Xj 3(k) a dot product is because the velocity vector is 
not directly observable for a specific spatial frequency. 
Even though the Kalman filter can generate estimates without 
the state being observable, estimation of the unobservable 
portion of the state relies on a priori information and not 
on the measurements  [Ref. 1]. Simplification of the 
individual filters results because only Xj;;(k) is estimated. 
The a priori information is later combined with all the 
estimates of X;;,(k) to yield a final velocity estimate. The 
final velocity estimate is calculated using the weighted least 


squares algorithm given by: 


$(k) = (w'z'w) w'z رم یں‎ (3.12) 


where w was defined by Eqn. 2.10 and 


1 
0 
P Ë 33 | 


1 NES) 


0 
هی ظ 
3,3 


9 , 


F IS the 3,3 component of the ith estimation error 


covariance matrix computed by each of the SFEKF’s. 
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Figure 3.1 The Modified Extended Kalman Filter 
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The state equation for the SFEKF is 


X,(K+1) = A,(v)X,(k) + Z,(K) 


D (v) O 
۱ 0 -1 X (k) + Z,(K) 3.14 


and Z, has the known covariance matrix 


6,۶ 0 0 
T 2 
Q2 2, ^ F[2 2  ] = 0 o; 0 ۰ (3.15) 


2 


The measurement equation and corresponding measurement 
noise covariance matrix for the SFEKF are similarly derived 
for each spatial frequency and have the form shown earlier in 
Eqn. 2.19. 


The linearized state equation is given by 
6X,(k*1) = À, (X, ) 8X, (K) + Z (k) , (3.16) 


where À, (X, )is the Jacobian for a specific spatial frequency 


and is linearized about X . 
۱ 


The extended Kalman filter equations, Eqns. 3.3 through 


3.8, are similarly modified to accommodate the SFEKF model. 


D. THE REDUCED ORDER EXTENDED KALMAN FILTER 

Data compression of images is achieved by truncating 
Spatial frequencies [Ref. 1]. The parallel structure of the 
MEKF can easily accommodate this truncation and suggests a 


reduced-order filter. Even though this will lead to a 
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deterioration of the output image quality, the degrading 
effects of the image noise are simultaneously reduced. 

The actual selection of, which spatial frequencies to 
truncate will depend upon the distribution of the energy in 
the signal and noise level. A priori knowledge of the spatial 
frequency content would assist in the selection of which 
Spatial frequencies are to be analyzed when filtering a series 


of image frames. 


18 


IV. SIMULATIONS 

The MATLAB program written for this thesis simulated an 
image frame with the dimensions of 32x32 pixels and is 
contained in Appendix C. The two types of moving objects 
employed for test purposes were an 8x8 pixel square of unity 
amplitude (Figure. 4.1) and an 8x8 pixel pyramid with 
amplitude intervals of 0.25 (Figure 4.2). Additionally, tests 
were conducted using a fixed a checkerboard background. 
Figures 4.3 and 4.4 depict the test objects moving across a 
checkerboard pattern. The individual checkers were 8x8 pixel 


squares with an amplitude of 0.25. 


A. PERIODICITY OF THE TEST OBJECTS 

An interesting challenge had to be overcome because of 
the small dimensions of the image frame. The circular shift 
property of the discrete Fourier transform allowed the 
computer simulation to be carried out over several hundred 
iterations. The original 32x32 pixel image frame used in the 
Simulations was assumed to be a periodic two-dimensional 
sequence. The moving object quickly traversed the image frame 
because of the small dimensions. By assuming periodicity, 
test objects that went off one edge of the frame would 
reappear on the opposite edge. This assumption allowed 


Simulations to be run over several hundred iterations because 
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tervals of 0.25. 


Figure 4. 


in 


An 8x8 pixel, pyramid object with amplitude 


2 





amplitude. 


square object of unity 





Figure 4.1 An 8x8 pixel, 


— 





Figure 4.3 An 8x8 pixel, square object of unity 
amplitude moving across a checkerboard background where 
each of the individual checkers has an amplitude of 

. .25. 
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Figure 4.4 An 8x8 pixel, pyramid object with amplitude 
intervals of 0.25 moving across a checkerboard 
background where the individual checkers have an 
amplitude of 0.25. 
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it appeared to the extended Kalman filter that the image frame 


was of a larger dimension. 


B. DISCRETE FOURIER TRANSFORM OF MOVING IMAGES 

A description of how the transformed image frame appears 
in the spatial frequency domain will help demonstrate why the 
extended Kalman filter was implemented in the spatial 
frequency domain. 

As an example, the square object depicted in Figure 4.1 
is transformed to the spatial frequency domain and the 
resulting magnitude plot is shown in Figure 4.5. Just as a 
one-dimensional square wave transforms to a one-dimensional 
Sinc function, the two-dimensional square object transforms 
to a two-dimensionai sinc function. 

What is not obvious is that as the square object is moved 
across the image frame, a linear phase shift results in each 
of the spatial frequencies; the magnitude plot remains the 
same. The amount of phase shift is frequency dependent as 
was shown in Equation 2.21. 

Figure 4.6 shows the linear phase shift that occurs in 
the spatial frequencies F(1,0) and F(0,1). When the square 
object is moved horizontally 1 pixel per iteration, a phase 
shift of -11.25° occurs in F(1,0). Likewise, when the square 
image is simultaneously shifted vertically 2 pixels per 
iteration, F(0,1) experiences a 22.5° phase shift per 


iteration. 
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Figure 4.5 Magnitude plot of the transformed image 
frame depicted in Fig. 4.1. 
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Figure 4.6 A linear phase shift of -11.25° per iteration 
occurs in F(1,0) (dashed line) and 22.5° per iteration for 
F(0,1) (solid line) when the square object is moved 1 
pixel horizontally and 2 pixels vertically per iteration. 
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Appendix B further discusses the idea of how objects 
moving at a constant velocity in a no-noise environment create 
a linear phase shift in the spatial frequency domain. Also 
shown is the effect of adding white noise to the image frame 


before it is transformed to the spatial frequency domain. 


C. LOWPASS FILTERING OF THE TRANSFORMED IMAGES 

After the test objects were transformed to the spatial 
frequency domain, the specific spatial frequencies were 
analyzed. Most of the energy for the two test objects was 
centered about the dc component which was located at the 
center of the transformed array. 

Because most of the energy was found near the dc 
component, the high frequency components were truncated or 


filtered out using a lowpass filter, 


Y(K,,k)) = X(k,,k,)H(k,,k,) (4.1) 


where Y(k,,k2) is the filtered image, X(k,,k>) is the 
transformed image, and H(k4,,k;) is the lowpass filter. 

As an example, Figure 4.1 was transformed to the spatial 
frequency domain where a lowpass filter was applied. This 
operation truncated the high frequency components. After 
transformation back to the spatial domain, Figure 4.7 shows 


the resulting filtered image frame. This filtered image frame 


Zo 
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Figure 4.7 The square image depicted in Figure 4.1 
after being filtered in the spatial frequency domain. 


The sharp edges are lost due to the loss of the high 
frequency components. 
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will be used in finding the error to calculate the Frobenius 


norm, which is discussed in the next section. 


D. PERFORMANCE MEASURES OF THE EXTENDED KALMAN FILTER 

After the specified number of iterations had been 
completed for each simulation, the final, filtered image frame 
was transformed back to the spatial domain. This filtered 
image frame could then be visually compared with the no-noise 
image frame to see how effective the extended Kalman filter 
performed. 

An additional performance measure involved analytically 
calculating the Frobenius norm throughout the simulation. 
The Frobenius norm is a scalar which is computed by squaring 
each of the terms in a matrix, summing each of the squares, 
and then taking the square root of the final sum (Eqn. 4.1) 


[Ref. 7]. 


مس سد 


N. N š 
| A |. = 2 2 | ajl (4.1) 
= J = 


Two Frobenius norms were calculated for each of the 
simulations. The first Frobenius norm was calculated from 
the error matrix that was found by subtracting the current 
estimate of the image frame from the current, no-noise image 
frame. Because the extended Kalman filter only analyzed 


specific spatial frequencies, the current no-noise image frame 
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was lowpass filtered so that only those spatial frequencies 
analyzed by the extended Kalman filter remained. As was shown 
earlier, Figure 4.1 is an example of a current, no-noise image 
while its lowpass filtered counterpart is depicted in Figure 
4.7: The current estimate of the image frame was then 
subtracted from the lowpass filtered image frame and the 


second Frobenius norm was calculated for this error matrix. 


E. SIMULATION RESULTS 

The first set of computer simulations had a square test 
object traversing the image frame. The extended Kalman filter 
was evaluated under a variety of test conditions which are 
summarized in Table 4.1. 

The two types of backgrounds were a homogeneous background 
with zero intensity and a checkerboard background. The 
velocity in Table 4.1 ils described in parentheses where the 
first number is the horizontal velocity and the second number 
is the vertical velocity. Finally, the standard deviation of 
the white noise describes the noise level that was introduced 
to the image frame during the successive iterations. 
Additionally, to ensure that the same random noise sequence 
was used for each of the conputer simulations, a random number 
seed was used to initilize the random number generator. 

Figures 4.8 through 4.47 are the computer simulation 
results using the square test object. The no-noise image 


frame, that the extended Kalman filter is estimating, is 
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frame that the extended Kalman filter is estimating is 
presented first. The extended Kalman filter estimate is then 
depicted; followed by the horizontal and vertical velocity 


estimates and the Frobenius norm plots. 


TABLE 4.1 SIMULATION TEST CONDITIONS (SQUARE OBJECT) 





H - homogeneous background, C - checkerboard background 

The results of the first two simulations are shown in 
Figures 4.8 through 4.15. The noise level for both 
simulations had a standard deviation of 1.0. The object was 
moved horizontally 1 pixel per iteration and simultaneously 
moved vertically 2 pixels per iteration. Even with the 
checkerboard background, the EKF was able to rapidly estimate 
the object velocity for both simulations. The Frobenius norm 
correspondingly reached steady state as the velocity estimates 


approached the actual object velocity. 
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The next two simulations are depicted in Figures 4.16 


607600 7ص‎ The standard deviation of the noise was 
increased to 2.0. Figures .4.16 and 4.20 show the square 
object position after 95 iterations. These two figures 


demonstrate the assumption that the image frame used in these 
simulations was a two-dimensional, periodic sequence. After 
95 lterations, the square object is leaving the image frame 
from the far right corner. “he sections of the square object 
that have left the image frame are shown reentering the image 
frame from the appropriate, opposite edges. 

The EKF was able to detect the square object motion as it 
traversed the homogeneous background. The high noise level 
caused the EKF to take considerably longer to estimate the 
object velocity, but the EKF was still able to provide an 
enhance estimate of the image frame (Figure 4.17). 

The increased noise level and addition of the checkerboard 
background precluded the EKF from detecting the square test 
object and estimating its velocity. 

The next set of simulations depicted in Figures 4.24 
through 4.31 were included to show how the EKF performed when 
the standard deviation of the noise was increased to 4.0. 
After 100 iterations, the EKF was unable to detect the square 
object or provide a reasonable estimate of the object 
velocity. 

The square object velocities of simulations #7 through 


#10 included noninteger values. The standard deviation of 
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the noise was kept at 1.0 for all four simulations. The 
response and performance of the EKF was similar to the first 
two simulations where the square object moved with integer 
velocities and the standard deviation of the noise was 1.0. 

The computer simulation results using the pyramid test 
object are depicted in Figures 40 through 71. The various 
test conditions for the pyramid object are summarized in Table 
m2. 

The pyramid test object contains slightly less than half 
of the signal energy content in the square test object. This 
had a direct effect on the performance of the EKF. 

Simulations #1 and #2 are shown in Figures 4.48 through 
4.55. The first two simulations using the pyramid test object 
had identical test conditions as the first two simulations 
using the square test object. However, the EKF needed about 
25 more iterations before it was able to track the motion of 
the pyramid test object (Figure 4.50). 

The addition of the checkerboard background in simulation 
#2 further reduced the signal energy of the pyramid test 
object. This adversely affected the performance of the EKF 
as it was unable to detect or track the motion of the pyramid 
test object. 

In simulations #3 and #4, the standard deviation of the 
noise was raised to 2.0. Here, the increased noise level 


prevented the EKF from detecting the pyramid test object. 
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TABLE 4.2 SIMULATION TEST CONDITIONS (PYRAMID OBJECT) 


(1.0,2.0) 












(1.0,2.0) 
)+) 2,0); 
(1.0,2.0) 
(2.5,1.5) 
(2.59155) 
(QOIS 


(1.0,0.5) 


H - homogeneous background, C - checkerboard background 


Simulations #5 through #8 had the pyramid test object 
moving at noninteger velocities. The standard deviation of 
the noise was kept at 1.0 for each of these simulations. The 
EKF was able to detect and closely estimate the noninteger 
velocities for simulations #5 and #7 which had the pyramid 
test object moving across the homogeneous background. Again, 
aS waS experienced in simulation #2, the addition of the 
checkerboard background precluded the EKF from detecting and 


tracking the pyramid test object. 
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for Simulation #1 after 40 Iterations. 
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Figure 4.9 Extended Kalman Filter Estimate of Image Frame 
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Figure 4.11 


Frobenius norm results for Simulation +77 
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Velocity estimates for Simulation #1. 
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Figure 4.12 Image Frame for Simulation #2 after 40 
Iterations. 
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Figure 4.13 Extended Kalman Filter Estimate of Current 
Image Frame for Simulation #2. 
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Frobenius norm results for Simulation #2. 
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Veclocity estimates for Simulation #2. 
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Figure 4.16 Image Frame for Simulation #3 after 95 
Iterations. 
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Figure 4.17 Extended Kalman Filter Estimate of Image 
Frame for Simulation #3 after 95 Iterations. 
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Frobenius norm results for Simulation #3. 
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Figure 4.18 Velocity estimates for Simulation #3. 
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Frobenius norm results for Simulation #4. 
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Velocity estimates for Simulation #4. 
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Image Frame for Simulation #5 after 100 
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Figure 4.25 


Frame for Simulation #5 after 100 Iterations. 
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Figure 4.26 Velocity Estimates for Simulation #5. 





Figure 4.27 Frobenius Norm Results for Simulation #5. 
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Figure 4.28 Image Frame for Simulation #6 after 100 


Iterations. 
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Figure 4.29 Extended Kalman Filter Estimate of Image 
Frame for Simulation #6 after 100 Iterations. 
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Figure 4.30 Velocity Estimates for Simulation £6. 





Figure 4.31 Frobenius Norm Results for Simulation #6. 
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Frame for Simulation #7 after 45 Iterations. 
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Frobenius Norm Results for Simulation #7. 
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Figure 4.34 


Velocity Estimates for Simulation #7. 
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Figure 4.36 Image Frame for Simulation #8 after 95 
Iterations. 
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Figure 4.37 Extended Kalman Filter Estimate of Image 
Frame for Simulation #8 after 45 Iterations. 
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Frobenius Norm Results for Simulation #8. 


Figure 4.39 
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Flgure 4.38 Velocity Estimates for Simulation #8. 
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Figure 4.40 Image Frame for Simulation #9 after 65 
Iterations. 





Figure 4.41 Extended Kalman Filter Estimate of Image 
Frame for Simulation #9 after 65 Iterations. 
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Figure 4.43 


Frobenius Norm Results for Simulation #9. 
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Figure 4.42 Velocity Estimates for Simulation £9. 
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Figure 4.44 Image Frame for Simulation #10 after 65 
Iterations. 
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Figure 4.45 Extended Kalman Filter Estimate of Image 
Frame for Simulation #10 after 65 Iterations. 
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Frobenius Norm Results for Simulation #10. 


Figure 4.47 
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Figure 4.46 Velocity Estimates for Simulation #10. 
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Figure 4.48 Image Frame for Simulation #1 after 40 
Iterations. 
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Flgure 4.49 Extended Kalman Filter Estimate of Image 
Frame for Simulation #1 after 40 Iterations. 
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Figure 4.50 Velocity Estimates for Simulation #1. 
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Figure 4.51 Frobenius Norm Results for Simulation #1. 
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Figure Es Image Frame for Simulation #2 after 40 
Iterati 
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Figure 4.53 xtende alm Filter Estimate of Image 
Frame for S "e tion E m er 40 ki erations 
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Frobenius Norm Results for Simulation #2. 
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Figure 4.54 


Velocity Estimates for Simulation #2. 
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Figure 4.56 Image Frame for Simulation #3 after 65 
Iterations. 
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Figure 4.57 Extended Kalman Filter Estimate of Image 
Frame for Simulation 43 after 65 Iterations. 
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Figure 4.58 Velocity Estimates for Simulation #3. 
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Figure 4.59 Frobenius Norm Results for Simulation #3. 
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Figure 4.60 Image Frame for Simulation #4 after 65 
Iterations. 
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Figure 4.61 Extended Kalman Filter Estimate of Image 
Frame for Simulation #4 after 65 Iterations. 
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Frobenius Norm Results for Simulation #4. 
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Figure 4.64 Image Frame for Simulation #5 after 95 
Iterations. 
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Figure 4.65 Extended Kalman Filter Estimate of Image 
Frame for Simulation #5 after 95 Iterations. 
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Figure 4.66 Velocity Estimates for Simulation #5. 
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Figure 4.67 Frobenius Norm Results for Simulation #5. 
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Figure 4.69 Extended Kalman Filter Estimate of Current 
Image Frame for Simulation #6. 
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Figure 4.71 Frobenius Norm Results for Simulation #6. 


Figure 4.70 Velocity Estimates for Simulation #6. 
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Figure 4.72 Image Frame for Simulation £7 after 65 
Iterations. 
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Figure 4.73 Extended Kalman Filter Estimate of Image 
Frame for Simulation #7 after 65 Iterations. 
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Figure 4.74 Velocity Estimates for Simulation #7. 
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Figure 4.75 Frobenius Norm Results for Simulation #7. 
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Figure 4.76 Image Frame for Simulation #8 after 65 
Iterations. 


>= = 
سح 





y N z p TY > 
` `. +4 0 Q^ 





Figure 4.77 Extended Kalman Filter Estimate of Image 
Frame for Simulation #8 after 65 Iterations. 
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Frobenius Norm Results for Simulation #8. 
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Figure 4.78 Velocity Estimates for Simulation #8. 


Iteration 


300 


250 


200 


— — — —— — 
= 


150 


== 
سے 
e‏ 
9 ې 


100 


سے ے 
اع — — — 


— — — = — 


50 





— — — — 
w ہے‎ 
— — 


0.8 
0.4 
-0.6 0 


e +‏ 
ې و 9 


531 217711521 2۸ 


V. CONCLUSIONS 

The extended Kalman filter (EKF), as applied to the moving 
object model, was evaluated under a variety of test 
conditions. These test conditions included various nolse 
levels, two types of backgrounds, and integer and noninteger 
velocities. Inclusion of the Frobenius norm provided an 
analytical measure of the EKF’s performance. The performance 
of the EKF was directly linked with its ability to detect and 
track the test object’s motion. 

For simulations using a homogeneous background, the EKF 
was unable to detect and track the square test object when 
the standard deviation of the noise was raised to 4.0. 
Similar results occurred when the pyramid test object was used 
and the standard deviation was increased to 2.0. 

The detection of the test object by the EKF is also 
related to the signal energy content of the test object in 
relation to the energy contained in the image frame 
background. The addition of the checkerboard background 
increased the energy content of the background. As a result, 
the simulations showed that the detection threshold of the EKF 
was reduced. The EKF was precluded from detecting and 
tracking the square test object’s motion across a checkerboard 
background when the standard deviation of the noise was raised 
above 2.0. Even with a standard deviation of 1.0, the EKF was 


unable to detect and track the pyramid test object. 
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These computer simulations have demonstrated that the EKF 
algorithm can successfully operate in a  high-noise 
environment. These results are encouraging and warrant 


continued research. 
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APPENDIX A. TUTORIAL ON TWO-DIMENSONAL SPATIAL FREQUENCIES 


A. PROGRAM DISCUSSION 

The program SPAPICT.M begins by placing an 8x8 square 
object of unity amplitude in the center of the frame (Fig. 
A.1). The image frame is then transformed to the frequency 
domain where each spatial frequency component is defined by 


the discrete Fourier transform [ Ref. 2 ]. 


F (K I) - Y ]هه درد رماع‎ -j2n( 3ye + کہ‎ | | ç (A.1) 
Agel meet l 2 
where 0 < kí < N,-1 and 0 < k; < N;-1. 

The transformed image frame is sent through an ideal 
lowpass filter (LPF), l.e., the high frequency components are 
truncated and set equal to zero. The filtered image frame is 
then transformed back to the spatial domain where the sharp 
edges are lost due to the high frequency components being 
filtered out (Fig. A.2). 

The image frame is transformed into a matrix that exhibits 
conjugate symmetry.  Conjugate symmetry results when a real 
image is transformed into the frequency domain. The number 
of computations for each spatial frequency and its complex 


conjugate can be reduced using Euler'/s identity, 


el’ + e@!J* = 2 cos(w) , fae 2) 


1 
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the Spatial Domain. 


Figure A.2 Filtered Frame after being Transformed Back to 





Depicted on a 32x32 Pixel Frame. 


Image of Unity Amplitude 


An 8x8 Square 


Figure A.1 





where w = KinitKkin; for the two-dimensional case. The 
contribution of each spatial frequency and its complex 


congugate was computed using 


F(k,,k,) + F*(k,,k,) 
27 
2|F(k,, kj) |cos| 5 (n,k,+n,k,) + <F(k,,k,) | , (A.3) 


where N = 32. 

As an example, the spatial frequency F(-1,-1) is the 
complex conjugate of the spatial frequency F(1,1). In the 
reconstruction of the filtered frame in the spatial frequency 
domain, the contribution of a specific spatial frequency and 
its complex conjugate were calculated using Eqn. A.3 and then 
summed with the previous frame. 

Figures A.3-A.8 show the successive summation of each 
spatial frequency and its corresponding complex conjugate in 
the reconstruction of the filtered frame in the spatial 
frequency domain. Each figure depicting a specific spatial 
frequency has been scaled and phase shifted based on the 
Fourier coefficients that were computed using the DFT. 

Figure A.3 depicts the summation of the dc component and 
F(0,1) and its complex conjugate. Figure A.4 depicts the 
contribution of spatial frequencies F(0,2) and F*(0,2) which 
are then summed with the frame shown in Fig. A.3 to produce 


ENG. A.5. 
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Figure 3 Summation of the dc component, F(0,1), 0 
F*(O,1). 
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Figure A.4 Spatial Frequency F(0,2) 
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Thls process ilis continued with the contributions from the 
other spatial frequencles being computed and summed with the 
frame shown in Fig. A.5. Figure A.6 shows the spatial 
frequency F(2,1) which has two cycles in the n,-direction and 
one cycle in the n;-direction. Figure A.7 shows the summation 
of all the contributions up to this point. The final frame 


shown in Fig. 8 is the same as the frame depicted in Fig. A.2. 
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Final summation of all spectral frequencies. 


Figure A.8 





B. PROGRAM LISTING 
SPAPIC.M is a stand-alone MATLAB program that 


reconstructs a filtered image in the frequency domain. 


kc با با پا‎ e e با پا با با پل‎ k k k k k k k k k k k k k k k k k k k k k k * k k k k k k k k k k k k k k k k k k k k پا‎ k * 
VARIABLE DEFINITION 


pict : Frame with square image as it would appear in 

spatial domain 

8x8 square image with unity amplitude 

Frequency vector listing the pertinent 
Spatial frequencies needed to 
reconstruct the filtered image in the 
Spatial frequency domain 


img 
freq 


BFtpic : Frame after transformation to the spatial 
frequency domain 

freqpict : Frame after lowpass filtering in spatial 
frequency domain 


pictifft  : Lowpass filtered Frame after transformation 
back to spatial domain 


o? o? o AO NO NO AO AO AO AO AO AO AO AO OM OP oe 


KOK OK O O OK CC oC oC o o o CC CC CC با‎ k k k k k k k k k k k k k kk k kk k 'k k k kk k k k k k k k k * 


oo 


CREATE THE PICTURE 


img < ones(8); 

pictl = zeros(32); 

۱1 ۲14 : 211,14۶ 211[( = img; 
[ n1 n2 ] = size(picti); 


% CREATE FREQUENCY VECTOR 


n = 0; 
for i = 0:3, 
if i == 0, 
for j = 0:3, 
n = n + 1; 
freq(n,:) = [ (2*pi*i)/n1 (2*pi*j)/n2 ]; 
end 
else 
for j = -3:3, 
n = n 1; 
red BD TO 2۰۳۹/1 2*pi*j)/n2 ]; 
end 
end 
end 
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$ PLOT THE FRAME IN THE SPATIAL DOMAIN 

mesh (pict1) 

title(’Frame -- Spatial Domain’) 

meta spaplot 

$ TAKE FOURIER TRANSFORM OF FRAME 

fftpictl = fftshift(fft2 (pict1l) 7 

% PASS THE TRANSFORMED FRAME THROUGH A LOWPASS FILTER 


trnpict - f£tpictl1([14:20]J 11٦ 
trnpictl = trnpict(:); 


$ TAKE INVERSE FOURIER TRANSFORM OF FILTERED FRAME 
freqpict - zeros(32); 

freqpict([14:20],[14:20])-» trnpict; 

pictifft = ifft2(fftshift(freqpict) ); 

% PLOT THE FILTERED FRAME IN THE SPATIAL DOMAIN 

mesh (pictifft) 

title(’Lowpass Filtered frame in the spatial domain’) 


meta spaplot 
pause(5) 


% DISPLAY AND SAVE THE MAGNITUDE AND PHASE OF THE 5 
TRANSFORMED PICTURE 


diary spadiary.mat 


disp(’Lowpass Filtered Frame -- Magnitude and Phase’) 
zl = ۲ abs(trnpict1'); angle(trnpicer mi 
diary off 


% CONSTRUCT THE FRAME IN THE FREQUENCY DOMAIN USING 
% FREQUENCIES THAT WERE NOT TRUNCATED 
zl = 21): 1:17 ل‎ 
zl = flipud(z1l’); 
pict - zeros(32); 
for n = و ول‎ 
u = freq(n,1); 
V — freqí(n,2); 
mag = zl(n,1);/ 
phase = zl(n,2); 
x و‎ 0177 
۷ = 0:31; 
y ] » meshdonm(x,y):; 


( XX,y 
Z 2 * mag * cos(u*xx * v*yy + phase); 
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pict = pict + z; 

subplot(211), mesh(z) 

title(’Constructing truncated image in frequency ’.. 
‘domain’ ) 

xlabel(['f(',num2str(u),',"',num2str(v),')']) 


subplot(212), mesh(pict) 
xlabel([{’Adding f(’,num2str(u),’,’,num2str(v),’) and’... 
fe emuUM2ZStr( WW), ,-,numZ2str(v),’) to the ’... 
^picture']); 
meta spaplot 
pause(3) 
clg 


subplot (111) 


end 


cig, 


$ GET A PLOT OF THE FINAL IMAGE CONSTRUCTED IN THE FREQUENCY 
$ DOMAIN 


mesh (pict) 
title(’Final image constructed in the frequency domain’) 
meta spaplot 
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APPENDIX B. PPHSE.M 


A. PROGRAM DISCUSSION 

PPHSE.M is a stand-alone program written for use with 
MATLAB to demonstrate the spatial frequency phase shift that 
results from moving an object across an image frame in the 
spatial domain. 

PPHSE.M sequentially moves an 8x8 unity square across a 
32x32 image frame in the spatial domain for 50 iterations. 
Two image frames are transformed to the spatial frequency 
domain. The first image frame consists of a homogeneous 
background with no added noise. The second image frame has 
Zero mean, Gaussian white noise. 

The real and imaginary Fourier coefficients and the phase 
in degrees are recorded and plotted for both the no-noise and 
noise-corrupted image frames for the two spatial frequencies 
F(1,0) and F(0,1). Only the phase is worth noting since the 
magnitude of the transformed image frame remains constant. 
In the spatial frequency domain, the effect of moving the 
Square object across the image frame is a phase shift in each 
of the spatial frequency components. 

Figure B.1 shows the phase shift in degrees for the spatial 
frequency F(1,0). The effect of moving the square object 
across the image frame 1 pixel per iteration is a 11.25" phase 


shift per iteration. The dashed line represents the noise- 
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corrupted image frame. The phase shift is no longer perfectly 
linear, but it still tracks the no-noise phase shift rather 
closely. 

Figure B.2 shows the phase shift in degrees for the 
spatial frequency F(0,1). The effect of simultaneously moving 
the square object 2 pixels per iteration vertically is a 22.5" 
phase shift per iteration. The phase shift for the noise 
corrupted image frame is represented by the dashed line and 
the nonlinear phase shift again generally follows the no-noise 
phase shift. 

It should be noted that the phase shift due to the object 
motion is a function of the spatial frequency. In this 
example, F(1,0) and F(0,1) were chosen because they represent 
the spatial frequencies that purely describe the horizontal 
and vertical motion of the moving object. 

Figures B.3-B.6 sequentially show the magnitude of the 
real and imaginary Fourier coefficients for the spatial 
frequencies F(1,0) and F(0,1). The dashed line again 
represents the coefficients of the noise-corrupted image 
frame. 

The real and imaginary Fourier coefficients for each of 
the spatial frequencies were two of the three the states that 
were estimated by the extended Kalman filter described in this 


thesis. 
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B. PROGRAM LISTING 


$ PPHSE.M is a stand-alone program that demonstrates 
$ the spatial frequency phase shift that results when 
$ a moving object traverses an image frame in the 

$ spatial domain 


clear $ clear out all variables in memory 
clear functions 


shpe - ['Square'/]; $ Define the shape of the moving image 

Vrow = 2.0; % Vertical Velocity -- # of pixels/its 

vclm - 1.0; $ Horizontal Velocity -- # of pixels/its 

imax - 3 $ Size of the spatial frequency truncation 

jmax = B 

its - 50; $ Number of iterations 

sigl = 1.0; % Standard deviation of the noise 

shape = ones(8); % Define the shape -- square of % 
unity amplitude 

pict = zeros(32); $ Create digital array to place % 

image in 


Establish normal distribution 
Define a standard random seed 
Define the size of the image 


rand('normal'); 
rand(’seed’ ,2380849) 
[ nl n2 [ = size(shape); 


oe o9 oe 


$ CREATE STORAGE VECTORS FOR CALULATED DATA 


phaseO-zeros(2,its); $ Phase shift in degrees -- no noise 
phasel-zeros(2,its); $ Phase shift in degrees -- 1.0 Std dev 
rsltO - zeros(2,its); % Fourier coef. -- no noise 

rsltl = zeros(2,its); % Fourier coef. -- Std dev. = 1.0 


% MOVE THE IMAGE ACROSS THE DIGITAL ARRAY FOR THE NUMBER OF 
% ITERATIONS SPECIFIED. 


for i - O:(its-1), 


$ DEFINE THE PARAMETERS FOR IMAGE LOCATION 


11 = rem (round (vrow*i+1),32) +1; 
12 = rem(round (vrow*i+8) ,32) +1; 
13 = rem(round(vclm*i+1),32)+1; 
14 = rem(round(vclm*i+8),32)+1; 
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$ PUT IMAGE INTO THE 32X32 ARRAY "pict" 
pict = zeros (32) ; 


1f 12 1 
11 4 5 55, 
pict(i1:32,13:32) 
۳1 2 2 7٦ 
۲1۷-۱ 0003 2 ۲٦ 
pict(1:12,13:32) 
else 
pict(il:32,1i3:14) 
pict(1:12,13: 114) 
end 
else 
if 14 > 
pict(11:12,13 ٦ 
pict(il:12, ۹۹۳ 
else 
۳ 1۳۲ (1 1: 1 م2‎ eres) 
end 
end 


shape(1:n1-i2,1:n2-i4); 

shape (n1+1-12:n1,n2+1-14:n2) ; 
shape(1:n1-i2,n2-1-i4:n2); 
shape (nl+1-i2:n1,1:n2-i4) ; 


shape(1:n1-i2,1:n2); 
shape(n1-*1-12:8,1:n2); 


shape(1:n1,1:n2-1i4); 
shape(1:n1,n2-41-i4:n2); 


shape(1:n1,1:n2); 


% CREATE THE FRAME CORRUPTED BY NOISE 
pictl = pict + (sigl*rand (32) ); 
% TRANSFORM THE FRAMES INTO THE SPATIAL FREQUENCY DOMAIN 


fftshift(fft2(pict); 
fftshift(fft2(pictl))s; 


fftpicto 
ffcpict] 


% OBTAIN THE PHASE IN DEGREES FROM SPATIAL FREQUENCIES 
%  F(0,1) & F(1%9) 


57.2958*angle(fftpictO(16,17)); 
57.2958*angle(fftpict1(16,17)); 
57.2958*angle(fftpictO(17,18)); 
57.2958*angle(fftpict1(217,18)); 


phaseO(1,i-*1) 
phasel1(1,1+1) 
phaseO(2,i-*1) 
phase1(2,i-*1) 


H H H Il 


% OBTAIN THE FOURIER COEFFICIENTS FOR SPATIAL FREQUENCIES 
% F(0,1) & F(1,0) 


fftpicto( 16 ip p 
fftpict1(16,17); 
۶ ۲ (1 17 ر‎ 
fftpictl (17 13 


rs1t0(1,i+1) 
rsl1t1(1,1i+1) 
rslto(2,1-1) 
rslt1(2,i+1) 


end 
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% PLOT THE PHASE SHIFT FOR SPATIAL FREQ F(0,1) IN DEGREES 


( > 1:15 ۶ 

plot(n,phase0(1,:), 'w',n,phase1(1,:), '—-—-w') 
title(’Phase shift for image velocity of (1,2) ’) 
xlabel(’Phase shift for F(0,1) in degrees’) 
ylabel (’Degrees’ ) 

meta pphse 

pause 


* PLOT THE PHASE SHIFT FOR SPATIAL FREQUENCY F(1,0) IN 
% DEGREES 


plot(n,phase0(2,:),’w’,n,phasel(2,:),’--w’) 
title(’Phase shift for image velocity of (1,2) ’) 
xlabel('Phase shift for F(1,0) in degrees") 
ylabel (’Degrees’” ) 

meta pphse 

pause 


% PLOT THE REAL FOURIER COEFF FOR SPATIAL FREQ F(0,1) 


peer, real (rsbtO(1,:)) ,’w’ pn, veal(rsltl1(1,:)),’--w’) 
title(’Real Fourier coefficients for F(0,1)’) 
Xxlabel(’Iteration’), ylabel(’Magnitude’ ) 

meta pphse 

pause 


* PLOT THE IMAGINARY FOURIER COEFF FOR SPATIAL FREQ F(0,1) 


plot(n,imag(rslt0(1,:)),’w’,n,imag(rslt1(1,:)),’--w’) 
title(’Imaginary Fourier coefficients for F(0,1)’) 
xlabel(’Iteration’), ylabel(’Magnitude’ ) 

meta pphse 

pause 


* PLOT THE REAL FOURIER COEFF FOR SPATIAL FREQ F(1,0) 


mmou(n,real(rslt0O(2,:)),’w’,n,real(rsltl(2,:)),’--w’) 
title(’Real Fourier coefficients for F(1,0)’) 
xlabel(’Iteration’), ylabel(’Magnitude’ ) 

meta pphse 

pause 


* PLOT THE IMAGINARY FOURIER COEFF FOR SPATIAL FREQ F(1,0) 


plot(n,imag(rsltO(2,:)),'w',n,imag(rslt1(2,:)),'--w/) 
title('Imaginary Fourier coefficients for F(1,0)^") 
xlabel(’Iteration’), ylabel(’Magnitude’ ) 

meta pphse 

pause 
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APPENDIX C.  PMEKF.M 


PROGRAM LISTING 


PMEKF.M is a stand-alone MATLAB program that implements the 


modified extended Kalman filter (MEKF) for moving objects. 


* de e de e k k او‎ k k k k k k k k k k او پر ار‎ k k k k k k k k k k k k k k k k k k k k k k k k k ار ار‎ K K K K K K K K KOK 
VARIABLE DEFINITIONS 


string defining shape of moving object 
vertical velocity in # of pixels per iteration 
horiz. velocity in # of pixels per iteraiton 
rectangular dimensions of spatial truncation 
window; located at center of transformed frame 
number of iterations 

Standard deviation of the added noise 

flag indicating user's intentions of whether 
or not to include checkerboard background 
amplitude of individual checkers 

for successive runs; number indicating 
specific run 

image frame containing checkerboard bkgrd; 
moving image included later 

vector listing two-dimensional frequencies to 
be analyzed 

length of the freq vector 

Frobenius norm of error matrix; subtract 
filtered image frame from original, no-noise 
image frame 
Frobenius norm 
filtered image 
image frame 
state estimate vector 

velocity estimate vector 

vector containing Fourler coefficients 

storage array for "pkk" covariance matrices 
estimation error covariance matrix; spatial 
frequency specific 

plant noise covariance matrix 

measurement noise covariance matrix 

no-noise image frame containing moving object 
no-noise image frame containing moving object 
and spatial frequency truncated 

image frame containing moving object and added 
zero-mean, white Gaussian noise 


of error matrix; subtract 
frame from truncated, no-noise 


shpe 
vrow 
vclm 


imax, jmax 


its 
sig 
chkrbd 


a 
runnum 


pictchkr 
freq 


lfreq 
normpic 


normpict 


xh 
vh 
Y 
pk 
pkk 


q 
E 


pict 
pictn 


randpict 


k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k با چاو از لر پر لا پار ما سلا‎ 


9 0 


oO AO oO oO AO oO oO oO oO oC 99 o€9 o9 o9 o9 o9 o9 o9 o9 o9 o9 o9 o9 o9 O9 o9 oO o9 o9 o9 oO o9 o9 o9 o? o9 o? o9 o9 o? 


ar 


cle 


clear functions 


shpe - ['Rectangle']; 
vrow = 2.0; 
vclm = 1.0; 
imax = 3 
jmax - 3; 
its = 2 ۶ 
slg = 1.00; 
chkrbd =  'n'; te wa -— yes:;eín' -- no ) 
a = 0.25; 
runnum = ; 
% DEFINE THE SHAPE 
xl = 1.00; x2 = 1.00; x3 = 1.00; x4 = 1.00; 
shape = [[ X1 x1 x1 x1 x1 x1 x1 x1 J] 
KIS K2 K2 X2 x2 x2 xl ] 
SNES X3 x3 x3 x2 xl ] 
[ X1 x2 x3 x4 x4 X3 x2 xl J 
[ xl x2 x3 x4 x4 x3 x2 xl ] 
0 1 فاد دص‎ X3 X3 X323 X2 X1 لا‎ 
lego x29x2 xX2 x2 x2 x2 x1-] 
[ xl xl xl xl xl x1 x1 x1 J); 
$ CREATE A CHECKERBOARD BACKGROUND IF REQUIRED 
if (chkrbd -- 'y'), 
pictchkr - zeros(32); 
chkr — a * ones(8); 
for i = 1:4, 

l f 1 == 1 == P 
pictchkr([i*8-7:i*8],([1:8]) » chkr; 
pictchkr((u*8-7:1*8),[17:24]) = chkr; 

else 
pictchkr([i*8-7:1*8],[9:16]) = chkr; 
pictchkr([i*8-7:i*8],[25:32]) » chkr; 

end 
end 
end 


E 


$ CREATE FREQUENCY VECTOR 


n = 0; 
for i - O:imax, 
if i == O, R 
for j = 0:jmax, 
n = n+1; 
freq(n,:) -E 
end 
else 
fer j = -jmax:jmax, 
n = + 
fred(n;:) m 1 ٣ 
end 
end 
end 


$ CREATE AND DEFINE NECESSARY VARIABLES 


SS 1/77 

jJ sqrt(-1); 

lfreq = length (freq) ; 
rand (‘normal’) ; 
rand(’seed’ ,2380849) 


Establish normal distribution 
Define a standard random seed 


o? o9 QO 


[ nl n2 ] = size(shape);% Dimensions of the image frame 
y = zeros(lfreq,its); 

xh = zeros(3,lfreq); 

vh = zeros(2,its); 


q= |] 4 O O 
0 0 O 
0 0.477 
r — sig^2*512*eye(2); 


normpic - zeros(1,its); 
normpict- zeros(1,its); 


$ GENERATE THE INITIAL ESTIMATION ERROR COVARIANCE MATRIX. 
$ EACH INDIVIDUAL FREQ COVARIANCE MATRIX IS STACKED IN "pk" 


pk = zeros(9,1freq) ; 
for k = i) freq. 


pkk = [ 1024 0 0 
0 1024 0 
0 0 4*freq(k,:)*freq(k,:)^']; 


pk(:,k)= pkk(:); 
end 
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$ MODIFIED EXTENDED KALMAN FILTER 


$ OBJECT IS MOVED ACROSS THE IMAGE FRAME FOR THE NUMBER OF 
% ITERATIONS SPECIFIED 


for i = 0: (its-1), 


% REINITIALIZE THE TYPE OF IMAGE FRAME USED IN THE 
% SIMULATION 


if (chkrbd == 'y'), 
pict = pictchkr; 
else 
pict 
end 


zeros(32); 


% DEFINE THE PARAMETERS FOR IMAGE LOCATION 


il = rem(round(vrow*i+1),32)+1; 
12 = rem(round(vrow*1+8),32)+1; 
i3 = rem(round(vclm*i+1),32)+1; 
14 = rem(round(vclm*i+8),32)+1; 


$ PUT IMAGE INTO THE 32X32 ARRAY "pict" 


1۶ 12 < il, 
ND 14 « i3, 
puct(11:32,.13:32) 
۱9 1:12, 1:14 ( 
۱۳ ۱۱ ۰ 2 1 14 } 
prct(1:3125323:32) 
else 
۱9 ۲ 11۰ م2‎ 1 3: 14( 
pict(1:12,1i3:14) 
end 
else 
٩1 4 > i3, 
pict(il:i2,i3:32) 
prct(il:i2,1:i4) 


shape (1:n1-12,1:n2-14); 

shape (n1+1-12:n1,n2+1-14:n2) ; 
shape (1:n1-12,n2+1-14:n2) ; 
shape (n1+1-12:n1,1:n2-i4) ; 


shape (1:n1-12,1:n2); 
shape (n1+1-12:8,1:n2); 


shape(1:n1,1:n2-i4); 
shape(1:n1,n2-1-i4:n2); 


I ll 


else 
pict(il:i2,i3:i4) = shape(1:n1,1:n2); 

end 
end 
randpict = pict + (sig*rand(32)); 
fpic = fft2(randpict); 
fpic = fftshift(fpic); 
fband = fpic((17-imax) : (17+imax) , (17-jmax) : (17+jmax) ); 
fvec سے‎ fband(:); 
y(:,i+1) = fvec(((length(fvec)+1)/2):length(fvec) ); 
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% TRUNCATE THE SPATIAL FREQS OF UNCORRUPTED PICTURE 


pictn سے‎ ۳ 7 

۳ 1 = fftshift(pictn); 

fband1 =pictn((17-imax):(17+imax),(17-jJjmax):(17+jmax)); 
pictn = zeros(32); 


1 


pictn((17-imax):(17+imax), (17-Jjmax):(17+j)jmax)) = fband1; 


pictn = fftshift(pictn); 
pictn — ifft2(pictn); 


$ ITERATE THRU THE "parallel" BANK OF EXT. KALMAN FILTERS 


for k = 1:17, 
dth = coxxh(3,K); 
cdth cos (dth); 
sdth sin(dth); 
dh = [ cdth sdth 
-sdth CdR TT 
dph = [ (-xh(1,k) *sdth+xh (2,k) *cdth) *co 
(-xh(1,k) *cdth-xh(2,k) *sdth) *co 


ah = [( dh dph 
0-0 TI, 
pkk(:) = pk(:,k)? 
pkplk = ah*pkk*ah’ +q; 
Gq = pkplk(: ,1:2)/ ) ۳۴۳۵) 1:2, OE 
xh(1:2,k) = dh*xh(1:2,K); 
yh xh(1:2 و(‎ 
ym = [ real(y Pr D) 
imag(y(K,i+1)) ] 
xh(:,k) = xh(:,k)+g*(ym-yh) ; 
pkk = pkplk - g* [1 0 0; O 1 Oj*pkplk; 
pk(:,k) = pkk(:); 


M 


end 


THESE THREE LINES GENERATE THE WEIGHTED LEAST SQUARES 
ESTIMATE OF THE VELOCITY VECTOR. COMMENTED OUT WHEN 
NOT USED. 


o? oo oe 


Sigin = diag(ones(1, (lfreq-1)) ./pk(9,2:lfreq)); 
fre = freq(2:lfreq,:)?; 
vh(:,i+1) = (fre’*sigin*fre) \fre’*sigin*xh(3,2:lfreq) ’; 


THIS LINE GENERATES THE LEAST SQUARES ESTIMATE OF THE 


% 

% VELOCITY VECTOR (unweighted). COMMENTED OUT WHEN NOT 
$ USED. 

% vh(:,i+1) = freq(2:lfreq,:)\xh(3,2:1lfreq)’; 
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۶ THIS LINE CHANGES THE ESTIMATES OF FREQUENCY FOR THE 
% INDIVIDUAL FILTERS LEAST SQUARES ESTIMATE. COMMENTED OUT 
% WHEN NOT USED. 


xh(3,:) = [{ freq*vh(:,i+1) ۰٠ 


% RECONSTRUCT THE IMAGE BACK INTO THE TIME DOMAIN. 


fpic = zeros(32); 

xhc - zeros(1,(2*lfreq-1)); 

xhc(lfreq:(lfreq*2-1)) -» xh(1,:)-jj*xh(2,:); 
xhc(lfreq:-1:1) - xh(1,:)-jj*xh(2,:); 

xhe(lfreq) = xh(1,1); 

fhand(:) = xhc; 
fpic((17-imax) : (17+imax) , (17-jmax) :(17+j}max)) = fband; 
fpic = fftshift(fpic) ; 

pic = ifft2(fpic); 


% CALCULATE THE NORM OF A MATRIX 

% FROBENIUS NORM : UNCORRUPTED IMAGE 

normplc(1+1) = norm(pict - pic,'fro')/ 32; 

$ FROBENIUS NORM : UNCORRUPTED, SPATIALLY-TRUNCATED IMAGE 

normpict(i+1)= norm(pictn - pic,'fro')/ 32; 

end 

$ DEFINE A TEXT LABEL FOR FINAL PLOTS 

name =[{’Std deviation = ',num2str(sig),' -- '... 
‘Iterations = ’,num2str(its),’ -- '... 
r'vel(',num2str(vrow),!','”",num2str(vclm),'!)'!]; 


% DISPLAY THE FINAL IMAGE 


pname = ['meta pmekf ',num2str(runnum) ]; 

clic 

mesh (pic) 

title([ shpe ’ -- Filtered Image -- Time domain']) 


xlabel(name) 
eval(pname) 
pause 
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mesh(pict) 

title([ shpe ^ -- Uncorrupted image’ }) 
xlabel(name) 

eval(pname) 

pause 


mesh (pictn) 

title([ shpe ’ -- Uncorrupted image (spatially truncated) ’}) 
xlabel (name) 

eval (pname) 

pause 


% PLOT FROBENIUS NORMS AND VELOCITY ESTIMATES 


n = 1:its; 
plot(n,normpic, ’-’,n,normpict, ’--’) 


title([ shpe ^ -- Frobenius Norm']) 
eval(pname) 

pause 
plot(n,vh(1,:),^9* nh vat. eee 
title([ shpe ’ -- Velocity Estimates’ ]) 


xlabel (name) 
eval (pname) 
pause 


% SAVE THE CALCULATED DATA 


m = ] n’ vh’ normpic’ normpict’]; 
ekfdata-['^save pmekf '/,num2str(runnum),^.mat m /ascii']; 
eval(ekfdata); 
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